arXiv: 1503.05521vl [cs.CV] 18 Mar 2015 


l 


Nonparametric Detection of Nonlinearly Mixed Pixels 
and Endmember Estimation in Hyperspectral Images 

Tales Imbiriba* 1 Student Member, IEEE, Jose Carlos Moreira Bermudez^, Senior Member, IEEE 
Cedric Richard^ 2 ), Senior Member, IEEE, Jean-Yves Tourneret ®, Senior Member, IEEE 

ll) Federal University of Santa Catarina 
Department of Electrical Engineering, 88040-900, Florianopolis, SC - Brazil 
tel.: +55.48.3721.7719 fax.: +55.48.3721.9280 

talesim @ gmail .com j .bermudez @ ieee. org 

Laboratoire Lagrange 

Universite de Nice Sophia-Antipolis, France 
phone: (33) 492 076 394 fax: (33) 492 076 321 

cedric.richard@unice.fr 

(3) Universite de Toulouse, IRIT/INP-ENSEEIHT/TSA, France 
phone: (33) 534 322 224 
jean-y ves .tourneret @ enseeiht.fr 


Abstract 

Mixing phenomena in hyperspectral images depend on a variety of factors such as the resolution of observation devices, the 
properties of materials, and how these materials interact with incident light in the scene. Different parametric and nonparametric 
models have been considered to address hyperspectral unmixing problems. The simplest one is the linear mixing model. Never¬ 
theless, it has been recognized that mixing phenomena can also be nonlinear. The corresponding nonlinear analysis techniques 
are necessarily more challenging and complex than those employed for linear unmixing. Within this context, it makes sense to 
detect the nonlinearly mixed pixels in an image prior to its analysis, and then employ the simplest possible unmixing technique to 
analyze each pixel. In this paper, we propose a technique for detecting nonlinearly mixed pixels. The detection approach is based 
on the comparison of the reconstruction errors using both a Gaussian process regression model and a linear regression model. 
The two errors are combined into a detection statistics for which a probability density function can be reasonably approximated. 
We also propose an iterative endmember extraction algorithm to be employed in combination with the detection algorithm. The 
proposed Detect-then-Unmix strategy, which consists of extracting endmembers, detecting nonlinearly mixed pixels and unmixing, 
is tested with synthetic and real images. 
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J.-Y. Tourneret was partly supported by ANR grants ANR-12- BS03-003 (Hypanema) and ANR-11-LABX-0040-CIMI. This work appeared in part in the 
Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP). Florence, Italy, May 2014 | lj. 
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I. Introduction 

Emerged in the 1960s with multispectral scanners, modern hyperspectral sensors produce two-dimensional hyperspectral 
images over a few tens to thousands of contiguous spectral bands El. Their high spectral resolution allows a comprehensive and 
quantitive analysis of materials in remotely observed data. This area has received considerable attention in the last decade. Due 
to historic downlink and computer processing limitations 0, hyperspectral images often trade spatial for spectral resolution 0. 
The observed reflectances then result from spectral mixtures of several pure material signatures. As a consequence, spectral 
unmixing has become an important issue for hyperspectral data processing 0. 

In a supervised setting, the spectral signatures of pure materials are available as vectors of reflectances of these materials at 
each wavelength. Such vectors are typically called endmembers due to their geometrical interpretation in the linear mixing case. 
Mixing phenomena depend on a variety of factors such as the resolution of observation devices, the properties of materials, 
and how these materials interact with incident light in the scene 0- Therefore, different parametric and nonparametric models 
have been considered to address hyperspectral unmixing problems. The simplest one is the linear mixing model, which assumes 
linear mixing of the endmembers contributions 0. It has been recognized that mixing phenomena can also be nonlinear 0,(6). 
The corresponding analysis techniques are necessarily more challenging and complex than those employed for linear unmixing. 
Nevertheless, nonlinear analysis of hyperspectral images has been widely explored in the past few years. See, for instance, 0, 
0-QE Nonlinear unmixing algorithms can lead to a better understanding of the individual spectral contributions, despite 
the increased complexity. Hence, it makes sense to detect the nonlinearly mixed pixels in an image prior to its analysis, and 
then employ the simplest possible unmixing technique to analyze each pixel. To this end, it is desirable to devise analysis 
techniques that combine endmember extraction, detection of nonlinearly mixed pixels and unmixing. 

The problems of extracting endmembers, detecting nonlinearly mixed pixels and unmixing are interlaced, and addressing 
them jointly is not a trivial task. For instance, most nonlinear unmixing techniques assume the endmembers to be known or 
to be estimated by an endmember extraction algorithm 11131 - 1121 II . However, most endmember extraction algorithms rely on 
the convex geometry associated with the linear mixing model (22 ) - |f26) , which obviously does not apply to nonlinearly mixed 
pixels. Endmember extraction techniques designed for situations where a significant part of the image is composed of nonlinear 
mixtures are rarely addressed in the literature. In fact, most of the techniques considering nonlinearly mixed pixels are part of 
a complete unsupervised unmixing strategy EE EE Detecting nonlinearly mixed pixels in an hyperspectral image is also 
a complex task. Physically motivated models £7), |[29l usually tend to be too complex for application in practical detection 
strategies. One possible approach is to consider a simplified parametric model for the nonlinearity. The parameters of this 
nonlinear model are then estimated from the image, and hypothesis tests are derived based on these estimates. For instance, 
a single-parameter polynomial post-nonlinear model is assumed in (30). The main question regarding parametric modeling of 
nonlinear mixing mechanisms is whether the chosen model can capture the actual nonlinear effects present in a scene. When 
nothing or little is known about the nonlinear mixing mechanism, a direct strategy is to exploit the property of linear mixing 
models to confine the noiseless data to a simplex. The hypothesis test proposed in lUl) is based on the distance between the 
observed pixel and this simplex. Though this test is robust to nonlinear mixing mechanisms, it conveys too little information 
about the nonlinearity as a tradeoff to guarantee simplicity. An alternative strategy is to use nonparametric techniques to 
extract information about the nonlinearity directly from the observations. A nonparametric unmixing technique based on kernel 
expansions is presented in (T3l , but this work does not address nonlinearity detection. A nonlinear mixing model for joint 
unmixing and nonlinearity detection is proposed in (32) . It assumes that the observed reflectances result from linear spectral 


3 


mixtures corrupted by a residual nonlinear component. This model is rather similar to the model initially introduced in Ifl3l , 
but the estimation method relies on a computationally intensive Bayesian procedure. 

All the detection methods discussed above assume known endmember spectral signatures. In most cases, the endmembers are 
assumed to have been estimated from the data. However, most endmember extraction algorithms exploit the convex geometry 
of the linear mixing model and assume the presence of pure endmember pixels in the image. They usually exploit one of the 
following properties: 1) the endmembers are the extreme points when projecting the data onto any subspace E2-E1, GI, 
EH, El; 2) the volume of a simplex spanned by any subset of points in the image is maximum when these points are the 
endmembers EH, El- Other strategies deal with nonnegative matrix factorization E3-E1- One method of particular interest 
to this work is the MVES (minimum volume enclosing simplex) algorithm l25l as it does not assume the presence of pure 
pixels, though it still exploits the geometry of linear mixtures. It solves a constrained least-squares optimization problem with 
a simplex volume regularizer. Practically, MVES finds the smallest simplex circumscribing the hyperspectral data. The vertices 
of this simplex are defined as the endmembers. 

In this paper, we propose a technique that combines endmember extraction and detection of nonlinearly mixed pixels in 
hyperspectral images. The detection approach is based on the comparison of the reconstruction errors using both a Gaussian 
process (GP) and a linear regression model. The two errors are combined into a detection statistics for which a probability 
density function can be reasonably approximated. We also propose an MVES-based iterative endmember extraction algorithm 
to be employed in combination with the detection algorithm to jointly detect nonlinearly mixed pixels and extract the image 
endmembers. The proposed method is tested with synthetic and real images. This work is organized as follows. Section [TT] 
reviews the linear mixing model and some nonlinear mixing models, and discusses different ways of modeling nonlinear 
interaction between light and endmembers. Section [III] discusses GP regression applied in the context of hyperspectral data 
unmixing. Section [TV] presents the application of GP to the detection of nonlinearly mixed pixels. Section [V] introduces a 
two-step iterative procedure to estimate the endmember matrix. This method combines the MVES algorithm and the nonlinear 
mixture detector proposed in Section [TV] Simulations with synthetic and real data are presented in Section [VT] Conclusions 
are finally presented in Section [VTI1 


II. Mixture models 

Each observed pixel can be written as a function of the endmembers plus an additive term associated with the measurement 
noise and the modeling error. Consider the model: 


r = ip(M) + n (1) 

where r = [n,..., Tl is a vector of reflectances observed in L spectral bands, M = [m i, ..., rn/p is the Lx R endmember 
matrix, whose *-th column rn, is an endmember, n ~ J\f(0, opI) is a white Gaussian noise (WGN) vector, and ip is an unknown 
mixing function. Several models of the form £[]) have been proposed in the literature, depending on the linearity or nonlinearity 
of ip, type of mixture, and other properties (6). 

A. The linear mixing model 

The linear mixing model assumes that each light ray interacts only with one material, disregarding multiple interactions 
between light and multiple materials 0 . The classical linear model assumes that ip is a convex combination of the endmembers. 
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In this situation, the vector r can then be written as 


r = Ma + n 

subject to 1 T a. = 1 and a >7 0 


( 2 ) 


where a = [a\,..., or] 1 is the vector of abundances of each endmember in M, R is the number of endmembers, and >; 
denotes the entrywise > operator. Therefore, the entries of a cannot be negative and should sum to one. The observation rp 
in the £-th wavelength of Q can be written as 

n = mj e a + n e (3) 


where m\ e denotes the f-th row of M as a column vector. In the noiseless case, namely, np = 0, the sum-to-one and positivity 
constraints over the abundances in © confine the data to a simplex. The vertices of this simplex are the endmembers, which 
justifies the terminology. 

Several parametric models have been proposed in the literature to describe nonlinear mixing mechanisms of endmembers 
in hyperspectral images. See ( 6 j and references therein. We shall now review two popular models that will be used later to 
generate synthetic data for evaluation purposes. 


B. Nonlinear mixing models 

The generalized bilinear model (GBM) 1T8I is given by 

fl-l R 

r = Mcx + Y J E 7 ij ctiCtj rrii © m 3 + n 

i=1 j=i+1 (4) 

subject to l T a = 1 and a >7 0 

where the parameters 7 y € [0,1] govern the amount of nonlinear contribution, and © denotes the Hadamard product. In the 
noiseless case, data following the model Q lie in a nonlinearly distorted simplex in whose vertices are the endmembers 
as in the linear case. For simplicity, here we consider a simplified version of this model in which the nonlinear contribution 
is controlled by a single parameter 7 such that 7 = 7 y for all (i, j). 

The post nonlinear mixing model (PNMM) {40l is given by 

r = g(Ma)+n (5) 

where g is a nonlinear function applied to the linear mixing model. The PNMM can represent a wide range of nonlinear 
mixing models, depending on the definition of g. For instance, the PNMM considered in fl~3) is given by 

r = (Ma)£ + n ( 6 ) 

where ( v)$ denotes the exponentiation applied to each entry of the vector v. For £ = 2, © becomes a bilinear model closely 
related to the GBM but without a linear term. The PNMM was explored in other works considering different forms for g 
applied to hyperspectral data unmixing S ED. 

The GBM Q and the PNMM © nonlinear mixing models mainly represent the scattering phenomenon where the light first 
interacts with an endmember, and then with a second one, before being captured by the hyperspectral sensor. Other models 
account for other kinds of interaction between light and endmembers, or consider other types of nonlinear effects. In the case 
of the intimate mixture model ll42l for instance, the endmembers are considered to be mixed at the molecular level. Other 
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nonlinear models can be considered depending on the characteristics of the scene (7|-|j9|, El, lH 8 l , ( 20 ), (29) , (40), (421 . 
More importantly, these informations are usually missing. Hence, it makes sense to develop nonparametric models that do not 
make strong assumptions about the type of nonlinearity involved in the mixture. 

III. Nonlinearity detection with Gaussian process regression models 
To detect nonlinearly mixed pixels in an hyperspectral image, assuming xp in (TJ is unknown, we propose to compare 
the reconstruction errors resulting from estimating xp with nonlinear and linear regression methods. Gaussian process (GP) 
regression methods consist of defining stochastic models for functions and performing inference in functional spaces l43l . The 
representation is rigorous but, at the same time, lets the data speak for themselves. This characteristic is desirable when little 
is known about the functions to be estimated. Using some knowledge obtained from the observations about the endmember 
matrix, we propose a supervised learning strategy to make inference on xp . 

This section describes the application of GP nonlinear regression to the problem at hand. Consider the training set {AT . r } 
with inputs AT = \rn Xl ..... rn, Xl ] ', and outputs or observations r = [r-|..... 77/ 1 . By analogy with the linear mixing 
model ( 0 ), we write the f-th row of {T} as 

re = xp ( m Xe ) + n e , (7) 


with Tf the £-th entry of the observation r, xp a real-valued function in a (reproducing kernel) Hilbert space TL, and rit an 
additive WGN in the f-th band. A Gaussian process is a collection of random variables, any finite number of which has a joint 
Gaussian distribution |43l . We define a Gaussian prior distribution for xp with mean and covariance functions given by 


E{xp(m Xe )} = 0 

E{xp(m Xt )xp(m Xer )} = n(m XlL1 Tn Xl ,) 


( 8 ) 


where re is a positive definite kernel. For notational simplicity, it is common but not necessary to consider GPs with a zero 
mean function. This assumption is not overly restricting as the mean of the posterior distribution is not confined to be zero 
(as shown by (ITTIl l. The prior on the noisy observation r becomes: 


r~AT(0,K + o*I L ), 


(9) 


with K the Gram matrix whose entries K tJ = K(m Xi ,m Xj ) are given by the kernel covariance function evaluated at m Xi 
and m Xj , up the noise power, and 7 the L x L identity matrix. 

To obtain the predictive distribution for 0* = xp(m Xt ) at any test point m Xt , we can write the joint distribution of the 
observation r and xp(m X:t ) as (43ll 


r 

ip * 

with re* = , m Xl ),..., rt(mx», mrj ] 1 and re** = The predictive distribution of -0*, or posterior of 

0*, is then obtained by conditioning (flOl i on the observation as follows: 


AM 0, 


K 


k: 


re* 

/‘V'jfcsI 


( 10 ) 


xp*\r,M,m x * ~ N (kJ [K + 1 r, re** - nj [K + ct^Il] 1 re*) . (11) 

The extension to a multivariate predictive distribution with test data AT* = \rn, Xr l ,..., rn, x * L ] 1 yields: 

xPJr, AT, AT* ~ N (Kj [K + a 2 n I L \ _1 r, T<T** - Kj [K + a 2 n I L \ ~ l TT*) (12) 
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where [if*]y = K(m Xiri , rn Xj ) and = K(m Xiri ,m Xtj ). Finally, we arrive at the minimum mean square error (MMSE) 

estimator for GP regression: 

V>* = E{V>*| r,M,M*} 

(13) 

= iC [A" + o^/i] r. 

In order to turn GP into a practical tool for processing hyperspectral data, it is essential to derive a method for estimating 
free parameters such as the noise variance <r 2 and possible kernel parameters defining the unknown parameter vector 0. We 
proceed as in f43l by maximizing the marginal likelihood p(r\M, <x 2 , 6) with respect to (a'p. 6), which leads to 

= argmax [K + a%I L \ 1 r - ^ log \K + a 2 n I L \^ • (14) 

This problem has to be addressed with numerical optimization methods. There is no guarantee that the cost function does not 
suffer from multiple local optima. However, our practical experience with hyperspectral data indicates that local optima are 
not a critical problem in this context. 

We conclude this section by introducing some kernel functions. Common examples include the linear kernel defined as: 


n{m Xi ,m Xj ) = mj.m x . 


(15) 


and radial basis function kernels, which depend on || m Xi — m Xj \\, such as the Gaussian kernel: 


K{m Xi ,m Xj ) = exp 



2 s 2 


\m X . 



(16) 


where s > 0 is the kernel bandwidth. In the sequel, we shall use the Gaussian kernel for its smoothness and non-informativeness, 
as we lack any knowledge about the unknown function ip. Note that this kernel has been used successfully in many signal and 
image processing applications, in particular for hyperspectral data unmixing E3, M- 


IV. Detection of nonlinearly mixed pixels 


A. The detection problem 

Given an observation r, we formulate the nonlinear mixture detector as the following binary hypothesis test problem 


Ho : r = Mol + n (17a) 

Hi : r = ip(M) + n (17b) 

where n is a zero-mean WGN with variance <x 2 . We assume that the endmember matrix M is available, or has been estimated 
from data using an endmember extraction technique J2J. We shall relax this hypothesis in Section [VJ and use the nonlinear 
mixture detector to jointly perform this task. 

We propose to compare the fitting errors resulting from estimating r with a linear or a nonlinear estimator ( fl3l ). Under 7f 0 , 
both estimators should provide good estimates. Under Hi, the estimation error resulting from the linear estimator should be 
significantly larger than that obtained with the nonlinear estimator. We shall now evaluate these fitting errors. 


B. Linear estimation error 

The MMSE estimator (fl3l > may be used with the linear kernel (fl5l) to estimate a. in (1 17 ab . Nevertheless, this would require 
to solve O in order to estimate ap. To save on unnecessary computational efforts, we shall limit the use of GP to nonlinear 
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model estimation. The MMSE estimator for ( 1 17 al> is given by: 

a = (M T M)~ 1 M T r (18) 

resulting in the following estimation error: 

eim = r - fii n = Pr (19) 

where P = 1^ — AT (' AT AT) “ 1 AT is an Lx L projection matrix of rank p = L — R. Note that no constraint is imposed on 
the abundance vector a. The objective is to obtain the best linear estimator, since the purpose at this point is not to perform 
unmixing, but to decide on the linearity (or not) of the considered mixture. 

Consider first the distribution for ||eij n || 2 . Under Hi, we have: 

e lia \Hi = P[i/> + n]. ( 20 ) 

This implies that 

e lin |Hi ~Af(Pt/>,a 2 P) (21) 

where we use that the projection matrix P is idempotent, that is, cr 2 PP T = a 2 P. Under Ho, we have: 

eii„|Ho ~Af(0,a 2 P). (22) 


Proper normalization of each squared entry ei; n ,i of en„ leads to the following conditional distributions under the two hypotheses: 

[Pjfp] 2 


c lin,i 


°nPi Pi 

P 2 

e lin,i 


vlpJPi 


Hi 


Xi 


< Pi Pi 


(23) 


Ho ~ x! (0) 


where pj denotes the ?-th row of matrix P, and % 2 (A) is the noncentral %-square distribution with n degrees of freedom 
and centrality parameter A Pffl . As P is idempotent and of rank p = L — R, which leads to ||eii n || 2 = r T Pr, we conclude 
that g5\ p. 33]: 


eiin 


Ho ~ x 2 P (0) • 


(24) 


C. Nonlinear estimation error with GP 

Since our interest at this point is not to make predictions for new data, but to evaluate the fitting error between the model 
output and the available data, we define the GP estimation error as: 


^nlin ^ ^ nlin 


where f n iin is given by (fl3l l with AT, = AT. Hence, using (fl3l) in (f25l) yields 


e n iin = 


= Hr 


M,=M 


1 -1 ■ 


where H = Il K t \K + cr?. I l\ ‘ is a real-valued matrix of rank L. 


(25) 


(26) 
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We shall now analyze the distribution of 11e n ii n 11 2 under hypotheses Ho and Hi. Under hypothesis Hi, we have: 


eniinl'Hi =H(ip + n). 


( 27 ) 


This leads to the following conditional distribution 

e n ii„|7fi ^U{H^,a 2 n HH T ). 

Under hypothesis Ho, the distribution for the error becomes 

e nhn \Ho ~ a 2 n HH T ). 

The distribution of the i-th entry of e n ] ln is thus given by 

eniin,i|~ J\f{hjMa , a 2 hjhi). 


(28) 


(29) 


(30) 


Proper normalization of each squared entry eniin,i of e n iin yields the following conditional distributions: 


e 1 

c nlin,i 

Hi r 


( [hjl>? \ 

cr 2 hj hi 

^cr 2 hj hi J 

e 2 

c nlin ,i 

n 0 r 

^xl 

([hjMcx ] 2 

cr 2 hj hi 

cr 2 hj hj 


(31) 


where hj denotes the z-th row of H. Non-central ^-square distributions in ( 1231 and (OH make the analysis of the test statistics 
in the next section intractable, even under Ho- In order to proceed, we argue that it is reasonable to assume that, under Ho, 
both the nonlinear GP regression method and the linear one should achieve the same level of accuracy. Considering d24l) . this 
approximation leads to 

n 0 = x 2 M (32) 

We validated this approximation using extensive Monte Carlo simulations. Figures [2ji and [2jt illustrate this assumption for a 
representative example. 


Klinll 2 

n-2 


D. The test statistics 


We propose to compare the squared norms of the two fitting error vectors e n ] ln and eu n to decide between Ho and Hi- Also, 
the test statistics should allow for the adjustment of the detection threshold to a given probability of false alarm (PFA) for 
design purposes. Considering these two objectives, we propose the following statistical test 


T = 


211 e n ii n | 


I e n lin II 2 + lleiinll 2 u 0 


(33) 


where r is the detection threshold. 

The reasoning behind the choice of T defined in (f33l) is as follows. Under Ho, both ||e n ii n || 2 and ||eii n || 2 are x-square 
dependent random variables. Now, we write e,\ m as e n i; n + \/2e, where e is assumed to be also a zero-mean i.i.d. Gaussian 
vecto|l and neglect the cross-term ej ]m e when compared to ||e || 2 in evaluating ||eii n || 2 under Ho- The latter approximation is 
due to the lack of correlation between e n | m and e, which can be largely attributed to mismatches resulting from the numerical 


1 The constant factor \/2 is for notation purpose only. 
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optimization required to solve (IT4b . Under these considerations, (l33l > can be written as T = ||e n iin|| 2 /(||e n iin || 2 + ||e|| 2 ) with 
both 11e n ii n 11 2 and ||e || 2 independent and x-square distributed. Such a statistics is known to follow a beta distribution ||46| . 

As the GP estimator tends to fit better nonlinearly mixed data, T should be less than 1 under hypothesis Hi- Conversely, T 
should be close to one for linearly mixed pixels, as ||e|| 2 tends to be much less than 2||e n ii n || 2 . Hence, as per (l33l >. we accept 
hypothesis Ho if T > r and we conclude for the nonlinear mixing hypothesis Hi if T < r. 

E. Determining the detection threshold 

Considering the assumption that the statistical test T has a beta distribution under Ho, a decision threshold r can be 
determined for a given PFA as 

r = H- 1 /J (PFA) (34) 

where B a j 3 is the cumulative distribution function of the beta distribution with parameters (a,/3). The parameters of this 
function must be estimated from the data. To this end, we initially determine an estimate A of the abundance matrix assuming 
the linear mixing model with the real observations R = [ri,..., rvy] and the known endmember matrix M. Then, using M 
and A we construct the synthetic image R s = MA, which satisfies Ho- For this linearly mixed hyperspectral image, we then 
compute, say, N samples of the test statistics T\Ho defined in (l33l > and fit a beta distribution to these samples. The threshold 
r for each PFA is then determined using ( l34l >. 

This procedure requires the knowledge of the endmember matrix M. The next section proposes an iterative technique to 
estimate M from an hyperspectral image, which we assume to contain linearly and nonlinearly mixed pixels. 

V. Endmember extraction in nonlinearly mixed hyperspectral images 

The presence of nonlinearly mixed pixels in an hyperspectral image tends to degrade the estimation accuracy of endmember 
extraction methods based on a linear mixing model. As a consequence, nonlinearly mixed pixels also affect the performance 
of algorithms using the endmember matrix such as the detection method presented in this paper. There has been few papers 
addressing endmember estimation from nonlinearly mixed images. A nonlinear unmixing algorithm is derived in l47l . The 
pixel reflectances are supposed to be post-nonlinear functions of unknown pure spectral components. A Bayesian strategy is 
proposed to both unmix the data and estimate the endmembers. Both tasks are however mutually dependent and the unmixing 
model is very specific. A nonlinear endmember estimation algorithm based on the approximation of geodesic distances is 
introduced in ll27l . 1481 . This algorithm can however suffer from the absence of pure pixels in the image, and the effectiveness 
of using manifold learning methods on real data still needs to be analyzed and confirmed. In this section, we propose an iterative 
technique for estimating the endmember matrix M under the reasonable assumptions that the number R of endmembers is 
known m-m, and that these endmembers are linearly mixed within at least a small part of the image. Nonlinear mixtures 
may however compose a significant part of the image. The proposed technique combines the detector of nonlinearly mixed 
pixels presented in Section [IV] and the endmember estimation algorithm MVES 125 1 . 

The procedure is described in Algorithm Q] It is a two-step iterative algorithm. The first step consists of using MVES to 
estimate the endmembers (line 2 and 14 in Algorithm Q3. The second step uses (l33l > to compute the detection statistics for 
all the L pixels in the image R lmp (line 7 in Algorithm H}. Then, all pixels whose detection statistic satisfies T(i) < T r are 
removed (line 9), where T r = rf x t (line 4 and 11) is the relaxed detection threshold. The use of a relaxed threshold is 
suggested to avoid discarding linear pixels during the first iterations, when the estimates of M are still not sufficiently accurate. 
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The relaxing factor is initialized for r/ = 0.9 and is increased by a factor r lnc = 0.1 /A' rnax at each iteration to improve pixel 
selection as the estimation of the matrix M improves (line 10). The procedure is repeated until the linear and the nonlinear 
GP models in (l33l) present similar fitting errors within the limit of e. Using this procedure, iy tends to the desired threshold r 
as the estimation of M improves, leaving mostly linear pixels for which both models have similar performance. A maximum 
number of iterations A' rnax is also set to avoid discarding too much data. 

Note that we have opted for the MVES algorithm for endmember extraction because it inscribes the data into a minimum- 
volume simplex. Thus, MVES is suitable to estimate M in the absence of pure pixels. This feature is specially interesting for 
our purpose since the procedure described above discards data, which may even be pure or near-pure pixels during the first 
iterations. Nevertheless, any other endmember estimation algorithm valid in absence of pure pixel could be potentially used 
with Algorithm Q] 


Algorithm 1: Iterative endmember estimation 


= 0.1/7V n 


PFA = 0.05; 


Input : The hyperspectral image R , and the number of endmembers R 
Output: Estimated endmember matrix M 

1 Initialization. T max — 1, T m i n — 0, s — 0.05, .Rtmp — R'i — 19? cc — 1, ry — 0.9, Ti nc 

2 M = MVES(-R tmp ,f?); 

3 Compute r using d34l >: 

4 

5 

6 

7 

8 
9 

10 
11 
12 
13 


T r = rf x r; %% (relaxed threshold) 

while T max Ti n i n > £ & cc <L N max do 
for 2 = 1 to Npixeis do 

Compute T(i ) using (133b : 

end 

Remove all pixels with T(i) < t v from f?, nlp ; 
r f = r f + ri nc \ %% (relaxing factor) 
r r = rf x t; 

T’max= rnax(T); T min = min(T); 
cc = cc+ 1; 


14 M = MVES(-R t m P , /?.); 

is end 


VI. Simulations 

This section presents simulation results to validate the proposed approach for detecting nonlinearly mixed pixels, with both 
synthetic and real images. The use of synthetic images is important as they provide a ground truth against which the performance 
of the detector can be verified. First, we propose a definition for a degree of nonlinearity of an hyperspectral image so that 
the relative performances of different detectors can be compared. This is helpful to quantify the relative energies associated 
with the linear and nonlinear mixing components in hyperspectral images generated with different nonlinear mixing models. 

A. Degree of nonlinearity 

Consider that a pixel vector can be written as the sum of a linear and a nonlinear mixing component as is the case for 
most existing nonlinear mixing models a, OH, ESI, ED: 

r = r hn + r nlin 

2 We do not account for noise contribution as it can be set by the user independently of the mixing model. 


( 35 ) 
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where r\ i n and r n ii n are, respectively, the linear and nonlinear mixing contributions to r. The energy of r is given by 

E = ||r || 2 = llr-^ll 2 + 2r 1 7 n r n]in + ||r nlin || 2 , (36) 


where En n = ||r'ii n || 2 is the energy of the linear contribution and E n ij n = 2rJ n r n n n + ||r n ii n || 2 is the part of the pixel energy 
affected by the nonlinear mixing. Given a mixing model, we define the degree of nonlinearity r]d as the ratio of the energy of 
the nonlinear contribution E n n n to the total energy E. Thus, 


Id 


/-Him 


1 

1 + A 


(37) 


where A = ||ri; n || 2 /(2rjTr n i; n + ||r n ]in|| 2 ). Next, we show how to apply this definition for generating synthetic samples with 
two different mixing models. 


1) Synthetic data generation with GBM: To be able to control the relative contributions of the linear and nonlinear mixing 
parts of the GBM model, we introduce a new scaling factor k into the generalized bilinear model (GBM) used in UTl . For an 
endmember matrix M and an abundance vector a, we write the modified noiseless GBM model as 


r = kM a + 


(38) 


where 0 < k < 1 , v = / '52f=i+i onotjm^Qmj is the nonlinear mixing term, 7 is the scaling parameter for the nonlinear 

contribution, and 0 is the Hadamard product. The degree of nonlinearity is then 

2k"f(L' T Ma) + 7 2 ||i'|| 2 

^ d k 2 \\Ma\\ 2 + 2ky(v T Ma) + 7 2 ||i / || 2 ^ 9 ) 

= 1 + A 

with A = k 2 \\Ma.\\ 2 /(2ky(v T Ma.) + 7 2 ||i'|| 2 ). We have to determine the scaling factors k and 7 so that the energy E is 
independent of rjd > 0 . This condition can be expressed as || Ma\\ 2 = fc 2 ||Ma || 2 + 2ky(y T Ma) + 7 2 ||i'|| 2 , leading to 


A = 


k 2 

1 - k 2 


(40) 


or _ 

k = \[rr^=V^Td. (4D 

To obtain 7, note that the denominator of A can be written as 7 2 |iz|| 2 + 2 k r ){u T Mcx.) = (1 — k 2 )\\Ma\\ 2 . Since 7 must be 
positive, we have 

7 = ( - 2 k{u T Ma) + ^ 4 k 2 (u T Ma ) 2 + 4||i/|| 2 (l - A; 2 )||Ma|| 2 ) . (42) 

Once k and 7 have been determined from r/d, we can generate data following the model in (l38l >. 


B. Synthetic data generation with PNMM 

To match the noiseless PNMM model © with the proposed formulation (l38l> . we complement it with a weighted linear 
mixture as follows: 


r = kMa + 7 u, 


( 43 ) 
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where v = (Ala.'f denotes the exponential value £ applied to each entry of Ala.. Model ( 143} reduces to ([ 6 ]) for k = 0 and 
7 = 1 . Again, parameters k and 7 are scaling factors that control the relative amounts of linear and nonlinear contributions 
given As for the GBM, both can be set using (TTTb and ( 142b . 


C. Simulations with known M 


We now present simulations with synthetic data and a known endmember matrix M. These simulations allow us to assess the 
detector performance disregarding estimation errors for the endmembers. Hence, they illustrate the potential of the proposed 
detector. To construct synthetic data, we used three materials (R = 3) extracted from the spectral library of the software 
ENVI™ |32): green grass, olive green paint and galvanized steel metal. Each endmember m r has L = 826 bands that were 
uniformly decimated by 3 to L = 275 bands. 

To evaluate the performance of the proposed detector, we generated 8000 synthetic samples by mixing the three collected 
spectra. Among the 8000 pixels, 4000 were generated using the linear model in ©■ and 4000 using the modified generalized 
bilinear model in (138} . A fixed abundance vector a = [0.6, 0.4, 0.1] T was used for all samples. Nonlinearly mixed samples 
were generated using different degrees of nonlinearity r/d £ {0.3, 0.5, 0.8} to test the detector under different conditions. The 
power of the additive Gaussian noise was set to <r^ = 0.001, which corresponds to SNR = 21dB. 

Figure U shows the receiver operating characteristics (ROCs) of the proposed GP detector and the LS robust detector presented 
in ED for the three values of rjd- The proposed detector performs better, especially for moderate to high degree of nonlinearity. 
For instance. Fig. [TJ: shows that the GP detector achieves a probability of detection of 1 for PFA = 0.1, while the FS robust 
detector yields a probability of detection of approximately 0.65 for the same PFA. Figure U shows the histograms of |e n | m || 2 , 

11 e ii n 11 2 and T for both linearly ('Hi,) and nonlinearly (Hi ) mixed data. The proposed test statistics clearly leads to histograms 
that differ significantly under both hypotheses Ho and Hi, which explains the improvement in detection performance. Figure [3] 
compares the histogram of T under Ho with the fitted beta distribution, confirming that the distribution of T can be reasonably 
approximated by a beta distribution. 

We considered two unmixing algorithms to assess the impact of the proposed detector on unmixing performance, one 
linear and one nonlinear. Finear unmixing was performed using the fully-constrained least-squares (FCFS) algorithm f53ll . For 
nonlinear unmixing, we used the SK-Hype algorithm II13 1 . The two algorithms were employed in two unmixing strategies. 
First, each algorithm was used to unmix the complete hyperspectral image. In the second strategy called Detect-then-Unmix, 
the proposed detector was used as a pre-processing step. Then, FCFS was used to unmix pixels detected as linearly mixed and 
SK-Hype was used to unmix pixels detected as nonlinearly mixed. The detection threshold r was determined for PFA = 0.01. 
Two synthetic images were considered with 1000 pixels each, 500 being linearly mixed and 500 being nonlinearly mixed. Each 
image was constructed using a particular nonlinear mixing model, with a fixed degree of nonlinearity r/d = 0.5 in both cases. 
The GBM (l38l > was used for the first image, while the PNMM {43} with £ = 3 was considered for the second image. The 
SNR was 21dB in both cases, and the abundances were drawn uniformly in the simplex. Parameters k and 7 were determined 
for each pixel vector to maintain the desired value of nonlinearity degree r)d for all simulations. To compare the results, we 
used the root mean square error (RMSE) of abundance estimation, defined as 


RMSE = 



N 


J2 H“'- 


2 


(44) 


where N is the number of pixels in each image. 
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(a) Robust LS detector ED- 
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(b) Proposed GP detector. 



(c) Comparison of the LS and GP detectors for r/d = 0.5. 

Fig. 1: Empirical ROCs for: (a) the Robust LS detector f3l! , (b) the proposed GP detector, (c) the two detectors for rjd = 0.5. 
All curves were obtained for 8000 pixels (4000 linearly mixed and 4000 nonlinearly mixed) and SNR = 21dB. Nonlinear 
mixtures were generated using the simplified GBM described in Section IVI-AI 
























14 


e n linl| 2 under Ho 



e nlinl| 2 under Hi 



ei in || 2 under Ho 



eu n || 2 under Hi 



T under Ho 



(a) Histogram of ||e 9 || 2 


(b) Histogram of ||ej|| 2 


(c) Histogram for the test statistics 03- 


Fig. 2: Histograms for (a) the squared norm of the GP fitting error, (b) the least-squares fitting error, and (c) the test statistics (l33l) . 


100 



| Test Statistics Histogram under H Q 
■ Adjusted Beta Distribution 
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Fig. 3: Histogram of the test statistics under 7fo and the adjusted Beta distribution. 


The results are presented in Tables [j] and |TT] For each image, these tables indicate the RMSEs for the linearly mixed part 
(LMM), for the nonlinearly mixed part (NLM), and for the full image (Full Img.) using the three unmixing strategies. The 
results in blue are those with the lowest RMSE in each row of the tables. As expected, FCLS has the best results when 
unmixing linearly mixed pixels. The same observation can be made for SK-Hype with nonlinearly mixed pixels. Nevertheless, 
we verify that the results using the Detect-then-Unmix strategy are very close to the best results for both types of pixels, LMM 
and NLM. When processing the whole image without prior information on the mixing nature of each pixel, the best results 
were those obtained with the Detect-then-Unmix strategy. 


D. Simulations with an unknown endmember matrix M 

The simulations conducted in Section lVI-Cl assumed the endmember matrix M to be known. Although this study is important 
to quantify the potential of the proposed detector, the endmembers are rarely known in practice. Hence, in this section, we 


TABLE I: RMSE in abundance estimation for M known and using the GBM mixing model (SNR = 21dB, r/ f / = 0.5). 



Image I: LMM + GBM 

Model 

FCLS 

SK-Hype 

Detect-then-Unmix 

LMM 

NLM 

Full Img. 

0.0173 ± 3.20e-04 

0.1220 ± 0.016268 

0.0871 ± 0.01362 

0.0349 ± 0.001543 

0.0644 ± 0.003463 

0.0518 ± 0.003054 

0.0189 ± 6.44e-04 

0.0666 ± 0.004054 

0.0490 ± 0.003548 
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TABLE II: RMSE in abundance estimation for M known and using the PNMM mixing model (SNR = 21dB, r} d = 0.5). 



Image II: LMM + PNMM 

Model 

FCLS 

SK-Hype 

Detect-then-Unmix 

LMM 

NLM 

Full Img. 

0.0170 ± 3.14e-04 

0.0637 ± 0.003696 

0.0466 ± 0.003228 

0.0354 ± 0.001723 

0.0551 ± 0.003745 

0.0463 ± 0.003047 

0.0173 ± 3.76e-04 

0.0559 ± 0.003636 

0.0414 ± 0.002946 


ROC Curve 



Fig. 4: ROCs for different proportions of nonlinearly mixed pixels and r)d = 0.5. Endmember extraction using VCA. 


study the sensitivity of the detection performance as a function of the endmember estimation accuracy and of the degree of 
nonlinearity. Endmember extraction was performed with the iterative method proposed in Section [V] and with VCA (23l for 
comparison. 

Figure0presents the results of 4 experiments using synthetic images with 5000 samples, SNR = 21dB, abundances uniformly 
sampled in the simplex, a proportion of nonlinearly mixed pixels in the image varying from 10% to 50%, and r)d = 0.5. For 
every experiment, the endmember matrix was extracted using VCA. These results show how the detection performance can 
degrade as the number of nonlinear pixels increases and as VCA loses accuracy in extracting the endmembers from the 
image. These results confirm the importance of devising alternatives to VCA (or tp other endmember extraction algorithms 
specifically designed for linearly-mixed images) for images containing nonlinearly-mixed pixels. Figure 0 presents the results 
obtained with Algorithm Q] for endmember extraction. For this experiment, we generated data with 50% of nonlinearly mixed 
pixels and different degrees of nonlinearity r/d £ {0.3,0.5, 0.8}. Comparing Fig. 0 and 0 shows that the results obtained with 
the iterative endmember extraction algorithm are very close to those obtained for a known endmember matrix M (which can 
be considered as the reference detector). 

Figurc0 illustrates a representative example of evolution obtained with the proposed iterative endmember extraction algorithm. 
These plots correspond to a simulation performed using 1000 synthetic samples, 500 being linearly mixed and 500 being 
nonlinearly mixed. The nonlinearly mixed pixels were created using the GBM (f38l) with r/d = 0.5. The data were projected 
onto the space spanned by the columns of the current endmember matrix M. They are represented as black dots. The current 
endmembers are shown as green dots. The true endmembers are shown as black circles at the vertices of the true simplex 
drawn with black lines. The data discarded at each iteration are shown within blue circles. Figure 0i shows the first iteration of 
Algorithm0 Numerous nonlinear samples are outside the simplex and endmember are poorly estimated. The situation improves 
in Fig. 03, which depicts the fourth iteration. Here, much less data lie outside the simplex, and two of the endmember estimates 
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ROC Curve 



Fig. 5: ROCs for different degrees of nonlinearity ij,i and 50% of nonlinearly mixed pixels in the image. Endmember extraction 
using Algorithm [I] 

have improved significantly. Similar improvement can be noticed in the seventh iteration in Fig. [(J. The final result obtained 
after 10 iterations only is shown in Fig. [6}l, where most of the nonlinear data were discarded and the endmember estimates 
are clearly close to the true endmembers. 

E. Real Data 

To test the proposed detector using real images, we used the data set available at the Indian Pines test site in North-western 
Indiana ED- This image was captured by the AVIRIS (Airborne Visible/Infrared Imaging Spectrometer). It has 145 x 145 
samples over 220 contiguous bands with wavelengths ranging from 366 to 2497 nm. Prior to analysis, noisy and water absorption 
bands were removed resulting in a total of 200 bands that were uniformly decimated to 50 to speed up simulations. The data 
set has a ground truth map that divides the samples into 16 mutually exclusive classes. In Table HIH the classes are organized 
by numbers (1 to 16), and the number of samples of each class is shown. Note, however, that the number of samples in 
each class can vary considerably. Note also that some classes are composed of different materials. We can count 20 different 
materials if we consider grass as an isolated material for the whole image. We chose to count each grass (depending on the 
accompanying material) as a different material, leading to 22 endmembers. Figures |7a| and [7b] display images from the Indian 
Pines region constructed by selecting three different bands, while Fig. [8a] presents the ground truth map for this image, where 
each class is represented by a different color. In Figure [8a] we also indice the class number for each area, where 0 represents 
the background, which is an unclassified area. 

To perform the simulations, we divided the image into eight sub-images to work with smaller areas of the image and to 
deal with 3 to 4 endmembers at a time. To define these sub-images, we also paid attention to balance the number of samples 
per endmember. By looking at Figs [7a] and [TbJ we can note that some classes seem to have materials that are not accounted 
for in the available ground-truth information. For instance, this is the case for classes 5,11 and 14. Therefore, we introduced 
extra endmembers for some of the sub-images. Table Hvl describes how the sub-images were organized, showing the classes, 
materials, numbers of pixels and endmembers chosen for each of the eight sub-images. 

For each sub-image, we estimated the endmembers as discussed in Section [V] with (V max = 10, a relaxing factor initially set 
to rf = 0.8, and incremented by r; nc = (1 — r/)/-ZV max = 0.2 at each of the 10 iterations. Then, we ran the detection algorithm 
with PFA = 0.001. We performed the unmixing step using FCLS for pixels detected as linearly mixed and SK-Hype for pixels 
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(a) First iteration. 


(b) Fourth iteration. 




(c) Seventh iteration. (d) Final result (after 10 iterations). 

Fig. 6: Graphical illustration of the endmember estimation process using the proposed iterative algorithm. The data set consists 
of 2000 pixels, with a proportion of 50% nonlinearly mixed pixels obtained with the GMB model and r]d = 0.5. Green dots 
are the current estimated endmembers, and black dots are the data projected onto the subspace spanned by the columns of the 
current matrix M. The true endmembers are shown as black circles at the vertices of the true simplex drawn with black lines. 
The data discarded at the corresponding iteration are shown within blue circles. 

detected as nonlinear mixtures. Figure [8b] presents the detection map superimposed to the ground-truth classes, where black 
dots represent pixels detected as nonlinearly mixed. 

Comparing the detection map in Fig. [8b] with Figs [7a] and [7b] one can note similarities between the detection map and some 
patterns observed in the image representations. For instance, the triangular shape in class 11 in Fig. [8b] is just besides what 
seems to be a road or trail when looking to Figure [7a] Similarities can be found between contours of detected nonlinear regions 
in Fig. [8b] and the corresponding regions in Figs [7a] or [Tb] Table [VI reports the RMSEs for the reconstruction error for each 
of the eight sub-images using three approaches, namely FCLS, SK-Hype, and Detect-then-Unmix. The results marked in blue 
correspond to the lowest RMSEs. For almost all sub-images, we note that the use of a nonlinear mixture detector improved 
the image reconstruction when compared to the pure linear or pure nonlinear unmixing strategies. 

VII. Conclusions 

This paper proposed a nonparametric method for detecting nonlinear mixtures in hyperspectral images. The performance of 
the detector was studied for supervised and unsupervised unmixing problems. Additionally, an iterative algorithm was derived 
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TABLE III: Indian Pines classes by region. 


Class number 

Class 

Num. of Samples 

1 

Alfalfa 

46 

2 

Corn-notill 

1428 

3 

Corn-mintill 

830 

4 

Corn 

237 

5 

Grass-pasture 

483 

6 

Grass-trees 

730 

7 

Grass-pasture-mowed 

28 

8 

Hay-windrowed 

478 

9 

Oats 

20 

10 

Soybean-notill 

972 

11 

Soybean-mintill 

2455 

12 

Soybean-clean 

593 

13 

Wheat 

205 

14 

Woods 

1265 

15 

Buildings-Grass-Trees-Drives 

386 

16 

Stone-Steel-Towers 

93 



20 40 60 80 100 120 140 20 40 60 80 100 120 140 

(a) Indian Pines representation (3-band combination #1) (b) Indian Pines representation (3-band combination #2) 

Fig. 7: Indian Pines test site representation selecting 3 different bands in (a), and 3 other bands in (b). 
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(a) Indian Pines ground truth. (b) Indian Pines detection map. 

Fig. 8: Detection of nonlinearly mixed pixels in Indian Pines hyperspectral image. Black pixels were detected as nonlinearly 
mixed ones by the proposed detector. 
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TABLE IV: Subimages organization 


Subimage 

Classes 

Materials 

# of pixels 

# of endmem. 

1 

9 and 7 

Oats and grass-pasture-mowed 

48 

3 

2 

1, 4 and 13 

Alfafa, wheat and corn 

488 

3 

3 

16 

Stone-steel-towers 

93 

3 

4 

15 

Buildings-grass-trees-drives 

386 

4 

5 

5 

Grass-Pasture 

483 

3 

6 

8 and 12 

Hay-windrowed and Soybean-clean 

1071 

3 

7 

3,6 and 10 

Corn-mintill, grass-trees and soybean-notill 

2532 

4 

8 

14 2 11 

Woods, corn-notill, soybean-mintill 

5148 

4 


TABLE V: Indian Pines recontruction error (RMSE) by subimage. 


Subimage 


RMSE ± STD 


FCLS 

SK-Hype 

Detect-then-Unmix 

1 

0.0028627 ± 6.6939e-06 

0.0030332 ± 6.0053e-06 

0.0029083 ± 6.5229e-06 

2 

0.0038963 ± 1.2293e-05 

0.003881 ± 9.4813e-06 

0.0038391 ± 1.1505e-05 

3 

0.0044259 ± 2.9087e-05 

0.0035981 ± 8.9722e-06 

0.0035537 ± 9.8622e-06 

4 

0.0040145 ± 1.1417e-05 

0.0039097 ± 8.0165e-06 

0.0038895 ± 8.5058e-06 

5 

0.0030848 ± 7.0516e-06 

0.0032353 ± 5.9761e-06 

0.0030527 ± 6.2275e-06 

6 

0.0039905 ± 6.5627e-06 

0.004055 ± 7.1531e-06 

0.0039644 ± 6.6603e-06 

7 

0.0034804 ± 5.8657e-06 

0.0035049 ± 5.9207e-06 

0.0034552 ± 5.9632e-06 

8 

0.0037665 ± 7.5723e-06 

0.0039314 ± 7.3092e-06 

0.0037531 ± 7.4932e-06 


for endmember estimation as a pre-processing step for unsupervised unmixing problems. It was shown that the combined use of 
the proposed detector and endmember estimation algorithm leads to better unmixing results when compared to state-of-the-art 
solutions. A degree of mixture nonlinearity based on the relative energies of the linear and nonlinear contributions to the 
mixing process was defined to quantify the importance of the linear and nonlinear model counterparts. Such a definition is 
important for a proper evaluation of the relative performances of different nonlinear mixture detection strategies. 
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